options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(1,4))
      drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
      drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
      drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
      drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

explanatory variable obs.x have noise

xN(x0.sx),yN(b0+b1*x0,s)

ex8-0.stan

normal regression without explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}

ex9.stan

normal regression with explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x10;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
  real<lower=0> sx;
  vector[N] x0;
}  
model{
  for(i in 1:N){
    x[i]~normal(x0[i],sx);
    y[i]~normal(b0+b1*x0[i],s);
  }
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x0[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] x1;
  vector[N1] y1;
  for(i in 1:N1){
    x1[i]=normal_rng(x10[i],sx);
    m1[i]=b0+b1*x10[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)

n1=10

#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -449.879 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       24       -33.767   0.000103114   0.000404158       0.825       0.825       28    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__     -33.77
##    b0         6.79
##    b1         1.89
##    s          3.28
##    m0[1]     15.48
##    m0[2]     27.72
##    m0[3]      9.45
##    m0[4]      2.88
##    m0[5]     29.60
##    m0[6]     37.52
##    m0[7]     19.87
##    m0[8]     34.79
##    m0[9]     11.39
##    m0[10]    31.28
##    m0[11]    29.52
##    m0[12]     9.94
##    m0[13]    26.71
##    m0[14]    46.07
##    m0[15]     8.25
##    m0[16]    21.68
##    m0[17]    25.69
##    m0[18]    18.13
##    m0[19]    41.38
##    m0[20]    34.53
##    y0[1]     16.95
##    y0[2]     29.28
##    y0[3]     15.41
##    y0[4]      1.50
##    y0[5]     29.19
##    y0[6]     38.80
##    y0[7]     14.03
##    y0[8]     32.00
##    y0[9]     19.81
##    y0[10]    32.30
##    y0[11]    30.46
##    y0[12]    17.10
##    y0[13]    29.38
##    y0[14]    46.54
##    y0[15]     4.22
##    y0[16]    24.14
##    y0[17]    25.40
##    y0[18]    16.10
##    y0[19]    45.04
##    y0[20]    36.64
##    m1[1]      2.88
##    m1[2]      7.68
##    m1[3]     12.48
##    m1[4]     17.27
##    m1[5]     22.07
##    m1[6]     26.87
##    m1[7]     31.67
##    m1[8]     36.47
##    m1[9]     41.27
##    m1[10]    46.07
##    y1[1]      3.97
##    y1[2]      5.35
##    y1[3]     14.07
##    y1[4]     15.00
##    y1[5]     25.49
##    y1[6]     28.47
##    y1[7]     33.38
##    y1[8]     33.03
##    y1[9]     39.83
##    y1[10]    44.76
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -34.23 -33.84 1.36 1.07 -36.89 -32.78 1.01      469      618
##    b0       6.92   7.01 1.57 1.40   4.25   9.60 1.01      211      215
##    b1       1.88   1.88 0.14 0.13   1.65   2.11 1.01      294      361
##    s        3.70   3.59 0.70 0.64   2.75   5.00 1.00     1158     1049
##    m0[1]   15.57  15.59 1.09 0.99  13.79  17.41 1.02      229      303
##    m0[2]   27.75  27.74 0.88 0.85  26.31  29.25 1.00     1297     1060
##    m0[3]    9.56   9.62 1.41 1.27   7.21  11.96 1.02      211      246
##    m0[4]    3.03   3.12 1.82 1.60   0.00   6.06 1.02      211      226
##    m0[5]   29.62  29.60 0.93 0.89  28.15  31.20 1.00     2021     1349
##    m0[6]   37.51  37.50 1.27 1.25  35.45  39.58 1.00     1789     1238
##    m0[7]   19.94  19.93 0.93 0.85  18.47  21.47 1.01      303      443
##    m0[8]   34.79  34.76 1.13 1.11  32.96  36.66 1.00     2220     1291
##    m0[9]   11.50  11.57 1.30 1.16   9.34  13.72 1.02      213      265
##    m0[10]  31.30  31.28 0.98 0.94  29.76  32.96 1.00     2337     1384
##    m0[11]  29.54  29.52 0.92 0.89  28.07  31.11 1.00     2030     1349
##    m0[12]  10.06  10.12 1.38 1.23   7.77  12.42 1.02      212      257
##    m0[13]  26.75  26.73 0.87 0.84  25.34  28.24 1.00     1023      992
##    m0[14]  46.01  45.99 1.78 1.72  43.08  48.92 1.00      818      981
##    m0[15]   8.37   8.45 1.48 1.32   5.89  10.89 1.01      211      218
##    m0[16]  21.75  21.73 0.88 0.84  20.37  23.22 1.01      381      517
##    m0[17]  25.74  25.71 0.86 0.82  24.34  27.18 1.00      811      960
##    m0[18]  18.21  18.22 0.98 0.92  16.65  19.86 1.01      260      350
##    m0[19]  41.35  41.33 1.49 1.46  38.92  43.75 1.00     1131     1077
##    m0[20]  34.53  34.50 1.12 1.10  32.73  36.40 1.00     2246     1291
##    y0[1]   15.48  15.38 3.95 3.81   9.23  22.04 1.00     1695     1585
##    y0[2]   27.91  27.80 3.89 3.58  21.65  34.44 1.00     1983     2009
##    y0[3]    9.62   9.62 4.12 3.92   3.05  16.36 1.00     1113     1322
##    y0[4]    3.18   3.20 4.25 3.88  -3.73  10.32 1.00      751     1086
##    y0[5]   29.61  29.60 3.85 3.76  23.36  35.95 1.00     1968     1927
##    y0[6]   37.48  37.59 3.95 3.69  30.97  44.04 1.00     2148     2057
##    y0[7]   19.92  19.96 3.87 3.69  13.56  25.93 1.00     1269     1702
##    y0[8]   34.60  34.64 3.89 3.69  27.88  40.94 1.00     1857     1933
##    y0[9]   11.44  11.45 4.03 3.79   4.80  18.13 1.00     1183     1681
##    y0[10]  31.27  31.30 3.97 3.89  24.83  37.70 1.00     1952     1943
##    y0[11]  29.55  29.58 3.81 3.45  23.22  35.77 1.00     2152     1933
##    y0[12]  10.03  10.00 4.02 3.87   3.44  16.48 1.00     1284     1763
##    y0[13]  26.64  26.67 3.88 3.80  20.16  32.95 1.00     1766     1851
##    y0[14]  45.98  46.07 4.16 4.10  39.10  52.82 1.00     1586     1655
##    y0[15]   8.36   8.39 4.02 3.78   1.75  14.81 1.00     1091      887
##    y0[16]  21.69  21.69 3.88 3.69  15.35  28.02 1.00     1749     1852
##    y0[17]  25.93  25.89 3.88 3.75  19.61  31.98 1.00     2043     1837
##    y0[18]  18.18  18.20 3.82 3.68  11.96  24.59 1.00     1763     1834
##    y0[19]  41.34  41.22 4.07 3.87  34.76  47.93 1.00     1970     1892
##    y0[20]  34.75  34.73 3.95 3.77  28.46  41.09 1.00     1973     1830
##    m1[1]    3.03   3.12 1.82 1.60   0.00   6.06 1.02      211      226
##    m1[2]    7.80   7.89 1.52 1.35   5.25  10.39 1.01      211      218
##    m1[3]   12.58  12.63 1.24 1.11  10.52  14.70 1.02      216      282
##    m1[4]   17.36  17.36 1.01 0.93  15.73  19.07 1.01      247      322
##    m1[5]   22.13  22.12 0.88 0.83  20.76  23.59 1.00      405      540
##    m1[6]   26.91  26.89 0.87 0.84  25.50  28.39 1.00     1063     1023
##    m1[7]   31.69  31.67 0.99 0.96  30.11  33.35 1.00     2355     1448
##    m1[8]   36.46  36.45 1.21 1.20  34.48  38.44 1.00     2006     1255
##    m1[9]   41.24  41.22 1.48 1.45  38.82  43.63 1.00     1141     1077
##    m1[10]  46.01  45.99 1.78 1.72  43.08  48.92 1.00      818      981
##    y1[1]    3.00   2.97 4.15 4.01  -3.72   9.75 1.00      908     1302
##    y1[2]    7.78   7.81 4.09 4.01   1.06  14.33 1.00     1116     1396
##    y1[3]   12.64  12.49 3.93 3.84   6.33  19.04 1.00     1433     1501
##    y1[4]   17.38  17.29 3.73 3.66  11.40  23.77 1.00     1746     1968
##    y1[5]   22.13  22.08 3.86 3.79  15.73  28.55 1.00     1928     1786
##    y1[6]   26.86  26.99 3.81 3.60  20.47  33.07 1.00     2065     1978
##    y1[7]   31.55  31.56 3.78 3.69  25.25  37.72 1.00     1964     2014
##    y1[8]   36.39  36.46 4.03 3.76  29.68  42.75 1.00     1988     1891
##    y1[9]   41.43  41.48 4.01 3.74  34.60  47.90 1.00     1851     1868
##    y1[10]  45.96  45.94 4.15 4.01  38.96  52.71 1.00     1899     1966

#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)

mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 1 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 2 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 3 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 4 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 4 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 1 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 2 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 3 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 4 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 4 finished in 0.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.6 seconds.
seeMCMC(mcmc,'[m,x]')
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -34.52 -36.49 11.42 10.37 -50.37 -12.77 1.06       73      122
##    b0       6.87   6.84  1.55  1.48   4.32   9.45 1.00      414      301
##    b1       1.88   1.87  0.14  0.14   1.64   2.13 1.01      402      276
##    s        2.70   2.81  1.15  1.21   0.77   4.51 1.03      120       38
##    sx       1.20   1.21  0.70  0.83   0.19   2.39 1.04      101      115
##    x0[1]    4.66   4.66  0.79  0.68   3.32   5.92 1.02      905     1262
##    x0[2]   10.62  10.74  0.91  0.79   8.99  12.00 1.01      618     1264
##    x0[3]    0.33   0.54  1.24  1.25  -1.97   1.89 1.01      287      381
##    x0[4]   -0.98  -1.22  1.22  1.27  -2.58   1.14 1.02      197      502
##    x0[5]   11.42  11.56  0.95  0.90   9.75  12.76 1.02      361      693
##    x0[6]   17.71  17.37  1.51  1.55  15.90  20.72 1.02      165       83
##    x0[7]    7.93   7.72  1.16  1.23   6.37  10.07 1.02      200      152
##    x0[8]   14.69  14.74  0.84  0.70  13.21  15.99 1.02      841     1998
##    x0[9]    1.83   1.99  0.97  0.82   0.05   3.11 1.01      433      453
##    x0[10]  14.20  13.94  1.32  1.38  12.54  16.67 1.02      169       61
##    x0[11]  12.09  12.07  0.83  0.69  10.75  13.49 1.01      582      740
##    x0[12]   2.49   2.30  1.06  1.05   1.05   4.35 1.02      221      416
##    x0[13]  10.01  10.15  0.91  0.84   8.42  11.35 1.01      541     1165
##    x0[14]  21.44  21.18  1.12  0.90  19.99  23.70 1.01      354      149
##    x0[15]   0.69   0.76  0.88  0.65  -0.86   2.05 1.01      990      598
##    x0[16]   7.29   7.44  0.92  0.85   5.65   8.59 1.01      415      615
##    x0[17]   9.50   9.64  0.89  0.82   7.95  10.80 1.01      511     1262
##    x0[18]   5.89   5.94  0.81  0.63   4.46   7.16 1.01     1511      813
##    x0[19]  18.43  18.35  0.92  0.71  16.99  20.15 1.01      518      223
##    x0[20]  13.56  13.75  1.18  1.30  11.49  15.09 1.02      228      505
##    m0[1]   15.63  15.66  1.54  1.39  13.14  18.17 1.00     1618     2603
##    m0[2]   26.78  26.81  1.76  1.75  23.99  29.59 1.01      523     2038
##    m0[3]    7.54   7.50  2.27  2.66   4.12  11.19 1.01      321     2267
##    m0[4]    5.05   5.18  2.50  2.76   0.80   8.66 1.02      183      283
##    m0[5]   28.27  28.34  1.84  2.02  25.34  31.14 1.02      295     1981
##    m0[6]   40.02  39.85  2.64  3.15  35.96  44.06 1.02      209     1656
##    m0[7]   21.73  21.60  2.15  2.48  18.41  25.10 1.02      218     1738
##    m0[8]   34.38  34.36  1.66  1.54  31.68  37.13 1.00     1679     2441
##    m0[9]   10.33  10.26  1.81  1.87   7.55  13.29 1.01      529     2619
##    m0[10]  33.45  33.31  2.34  2.76  29.93  37.09 1.02      232     1933
##    m0[11]  29.53  29.56  1.54  1.32  27.00  32.01 1.01     4331     2220
##    m0[12]  11.55  11.60  2.08  2.29   8.14  14.76 1.02      214      301
##    m0[13]  25.64  25.69  1.75  1.79  22.86  28.45 1.01      484     1835
##    m0[14]  47.00  47.15  2.00  1.92  43.59  50.08 1.00      985     2790
##    m0[15]   8.20   8.15  1.68  1.47   5.52  11.03 1.01     1428      797
##    m0[16]  20.55  20.61  1.74  1.81  17.69  23.32 1.02      355     1589
##    m0[17]  24.69  24.69  1.71  1.67  21.93  27.44 1.01      538     2261
##    m0[18]  17.92  17.91  1.53  1.37  15.48  20.45 1.00     1885     2411
##    m0[19]  41.38  41.41  1.75  1.60  38.48  44.19 1.00     2872     2521
##    m0[20]  32.28  32.41  2.34  2.78  28.65  35.92 1.02      209     1863
##    y0[1]   15.54  15.57  3.25  2.78  10.22  20.88 1.01     3309     3305
##    y0[2]   26.74  26.46  3.43  2.83  21.54  32.71 1.00     2157     3056
##    y0[3]    7.58   7.05  3.69  3.42   2.30  14.27 1.01      789     2073
##    y0[4]    5.07   5.60  3.85  3.60  -1.74  10.55 1.00      525     1308
##    y0[5]   28.32  27.99  3.53  3.16  23.11  34.50 1.00     1031     2844
##    y0[6]   39.99  40.48  3.95  3.98  32.78  45.51 1.01      477     3051
##    y0[7]   21.71  22.14  3.67  3.43  15.16  26.97 1.01      611     2353
##    y0[8]   34.35  34.20  3.32  2.85  29.06  40.02 1.01     3863     2907
##    y0[9]   10.37   9.95  3.46  3.01   5.15  16.45 1.01     1859     2480
##    y0[10]  33.47  34.01  3.76  3.54  26.82  38.68 1.01      557     2713
##    y0[11]  29.55  29.55  3.28  2.75  24.19  35.05 1.01     4263     3552
##    y0[12]  11.50  11.92  3.58  3.31   5.20  16.58 1.01      674     2425
##    y0[13]  25.61  25.24  3.45  2.98  20.36  31.68 1.00     2123     3397
##    y0[14]  46.92  47.25  3.60  3.09  40.74  52.38 1.00     2417     2529
##    y0[15]   8.21   8.10  3.45  2.83   2.61  14.05 1.01     3410     2279
##    y0[16]  20.52  20.21  3.42  3.05  15.33  26.39 1.00     1461     2592
##    y0[17]  24.75  24.41  3.33  2.86  19.81  30.67 1.00     1712     2853
##    y0[18]  17.94  17.80  3.33  2.82  12.48  23.55 1.01     3751     3456
##    y0[19]  41.35  41.42  3.45  2.79  35.59  47.10 1.00     3917     2795
##    y0[20]  32.24  31.67  3.78  3.48  26.98  39.13 1.01      558     1812
##    m1[1]    2.98   2.97  1.80  1.73  -0.01   5.95 1.01      409      308
##    m1[2]    7.75   7.71  1.50  1.43   5.29  10.24 1.00      414      302
##    m1[3]   12.51  12.50  1.23  1.22  10.52  14.58 1.00      417      325
##    m1[4]   17.28  17.29  1.02  1.00  15.69  19.01 1.01      421      494
##    m1[5]   22.04  22.06  0.91  0.86  20.51  23.52 1.02      412      115
##    m1[6]   26.81  26.83  0.94  0.88  25.08  28.35 1.02      447       56
##    m1[7]   31.57  31.61  1.09  1.04  29.62  33.34 1.01      418       60
##    m1[8]   36.34  36.36  1.33  1.28  33.96  38.54 1.01      390       86
##    m1[9]   41.10  41.09  1.62  1.58  38.37  43.82 1.01      386      186
##    m1[10]  45.87  45.85  1.93  1.89  42.74  49.12 1.01      385      450
##    x1[1]   -2.10  -2.09  1.36  0.90  -4.41   0.16 1.02     4210     2527
##    x1[2]    0.51   0.50  1.34  0.92  -1.71   2.70 1.02     4279     1693
##    x1[3]    2.98   3.00  1.39  0.92   0.69   5.29 1.02     3679     1846
##    x1[4]    5.56   5.55  1.38  0.90   3.32   7.87 1.02     3667     2424
##    x1[5]    8.09   8.09  1.40  0.92   5.78  10.47 1.02     4233     2071
##    x1[6]   10.62  10.63  1.36  0.88   8.34  12.79 1.02     3558     1978
##    x1[7]   13.14  13.18  1.39  0.93  10.76  15.38 1.02     3996     2192
##    x1[8]   15.72  15.71  1.39  0.94  13.49  18.08 1.02     4210     1766
##    x1[9]   18.26  18.23  1.38  0.92  16.07  20.66 1.02     3677     2009
##    x1[10]  20.76  20.79  1.38  0.90  18.36  23.01 1.02     3681     1907
##    y1[1]    3.08   2.93  3.53  3.15  -2.51   8.91 1.00     1225     1685
##    y1[2]    7.75   7.64  3.26  2.82   2.28  13.11 1.00     1422     2882
##    y1[3]   12.50  12.40  3.18  2.74   7.29  17.75 1.00     1989     2478
##    y1[4]   17.30  17.27  3.15  2.62  12.22  22.47 1.01     3246     1710
##    y1[5]   22.02  22.00  3.05  2.62  16.99  27.12 1.00     3624     3149
##    y1[6]   26.75  26.72  3.14  2.67  21.68  31.92 1.01     2488     2889
##    y1[7]   31.58  31.55  3.17  2.71  26.30  36.74 1.00     2193     2890
##    y1[8]   36.37  36.33  3.21  2.91  31.33  41.69 1.00     2054     3064
##    y1[9]   41.17  41.09  3.29  3.05  36.04  46.66 1.00     1656     2741
##    y1[10]  45.93  45.87  3.61  3.29  40.08  51.78 1.00      999     2168

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

outlier

ex10.stan

objective variable have outlier, y~cauchy(b0+b1*x,s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~cauchy(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=cauchy_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=cauchy_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)

n1=10

x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -1331.69 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       16      -32.5887    0.00192501   0.000150369           1           1       22    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -32.59
##    b0         7.16
##    b1         1.63
##    s          3.09
##    m0[1]     12.05
##    m0[2]     17.06
##    m0[3]     10.23
##    m0[4]     14.16
##    m0[5]      9.01
##    m0[6]     18.58
##    m0[7]      8.10
##    m0[8]     13.12
##    m0[9]     17.37
##    m0[10]    14.72
##    m0[11]     8.89
##    m0[12]    16.11
##    m0[13]    14.70
##    m0[14]    20.16
##    m0[15]    15.42
##    m0[16]    17.59
##    m0[17]    14.98
##    m0[18]    17.09
##    m0[19]    13.41
##    m0[20]    12.08
##    y0[1]     17.72
##    y0[2]     16.53
##    y0[3]      9.67
##    y0[4]     14.59
##    y0[5]     10.35
##    y0[6]     18.42
##    y0[7]      7.84
##    y0[8]     17.06
##    y0[9]     17.46
##    y0[10]    18.80
##    y0[11]     5.35
##    y0[12]    17.40
##    y0[13]    12.22
##    y0[14]    15.22
##    y0[15]    14.65
##    y0[16]    16.45
##    y0[17]     9.87
##    y0[18]    18.11
##    y0[19]    13.72
##    y0[20]    11.10
##    m1[1]      8.10
##    m1[2]      9.44
##    m1[3]     10.78
##    m1[4]     12.12
##    m1[5]     13.46
##    m1[6]     14.80
##    m1[7]     16.14
##    m1[8]     17.48
##    m1[9]     18.82
##    m1[10]    20.16
##    y1[1]      5.52
##    y1[2]     10.88
##    y1[3]     12.94
##    y1[4]     18.52
##    y1[5]     15.21
##    y1[6]     15.75
##    y1[7]     10.59
##    y1[8]     23.25
##    y1[9]     20.31
##    y1[10]    22.45
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -33.06 -32.71 1.33 1.12 -35.66 -31.62 1.02      384      618
##    b0       7.22   7.20 1.86 1.73   4.21  10.35 1.02      251      135
##    b1       1.61   1.63 0.39 0.37   0.95   2.23 1.01      309      385
##    s        3.50   3.42 0.65 0.60   2.61   4.66 1.00     1243     1508
##    m0[1]   12.07  12.08 0.96 0.91  10.46  13.66 1.02      325      306
##    m0[2]   17.01  17.01 1.04 0.98  15.26  18.73 1.00     1831     1364
##    m0[3]   10.27  10.27 1.24 1.15   8.23  12.39 1.02      262      137
##    m0[4]   14.15  14.16 0.80 0.75  12.80  15.46 1.00      871      887
##    m0[5]    9.06   9.05 1.48 1.36   6.67  11.60 1.02      252      109
##    m0[6]   18.52  18.54 1.30 1.26  16.34  20.60 1.00     1004     1270
##    m0[7]    8.15   8.14 1.66 1.53   5.46  10.94 1.02      250      121
##    m0[8]   13.12  13.14 0.85 0.79  11.68  14.49 1.01      456      443
##    m0[9]   17.32  17.33 1.09 1.03  15.50  19.07 1.00     1676     1345
##    m0[10]  14.70  14.70 0.81 0.74  13.35  16.04 1.00     1342     1026
##    m0[11]   8.93   8.92 1.50 1.38   6.51  11.51 1.02      252      107
##    m0[12]  16.08  16.09 0.91 0.84  14.54  17.59 1.00     2234     1403
##    m0[13]  14.69  14.69 0.81 0.74  13.34  16.02 1.00     1325     1026
##    m0[14]  20.09  20.09 1.61 1.60  17.42  22.62 1.00      697      967
##    m0[15]  15.39  15.38 0.85 0.78  13.99  16.79 1.00     2124     1323
##    m0[16]  17.54  17.56 1.12 1.07  15.64  19.34 1.00     1504     1343
##    m0[17]  14.97  14.96 0.82 0.77  13.60  16.31 1.00     1676     1151
##    m0[18]  17.05  17.05 1.04 0.99  15.28  18.77 1.00     1814     1364
##    m0[19]  13.41  13.43 0.83 0.78  11.98  14.75 1.01      528      577
##    m0[20]  12.10  12.11 0.95 0.90  10.50  13.68 1.02      327      306
##    y0[1]   11.96  11.99 3.60 3.48   5.90  17.92 1.00     1107     1363
##    y0[2]   16.92  16.93 3.69 3.65  11.05  22.87 1.00     1711     1926
##    y0[3]   10.35  10.42 3.69 3.64   4.35  16.44 1.00     1502     1802
##    y0[4]   14.17  14.18 3.71 3.53   8.06  20.14 1.00     1886     2011
##    y0[5]    9.18   9.16 3.86 3.69   3.06  15.42 1.00     1039     1682
##    y0[6]   18.55  18.67 3.80 3.64  12.33  24.85 1.00     1895     2014
##    y0[7]    8.12   8.18 4.01 3.85   1.53  14.90 1.00      996     1239
##    y0[8]   12.98  13.01 3.68 3.56   7.08  18.97 1.00     1768     1915
##    y0[9]   17.35  17.42 3.86 3.62  11.24  23.67 1.00     1656     1955
##    y0[10]  14.50  14.48 3.63 3.47   8.58  20.57 1.00     2033     1961
##    y0[11]   8.97   8.98 3.75 3.64   2.70  15.03 1.00      981     1730
##    y0[12]  16.25  16.21 3.71 3.42  10.19  22.38 1.00     2224     2008
##    y0[13]  14.79  14.77 3.71 3.54   8.63  20.84 1.00     1917     2055
##    y0[14]  20.16  20.08 3.99 3.82  13.76  26.88 1.00     1485     1699
##    y0[15]  15.35  15.33 3.64 3.52   9.28  21.06 1.00     1972     1839
##    y0[16]  17.52  17.50 3.91 3.75  11.20  23.88 1.00     1943     1893
##    y0[17]  14.79  14.75 3.67 3.51   8.70  20.82 1.00     1827     1890
##    y0[18]  17.08  17.03 3.83 3.71  10.88  23.25 1.00     2067     1971
##    y0[19]  13.40  13.35 3.70 3.63   7.31  19.80 1.00     1974     1904
##    y0[20]  12.00  12.00 3.72 3.51   5.96  17.98 1.00     1526     1745
##    m1[1]    8.15   8.14 1.66 1.53   5.46  10.94 1.02      250      121
##    m1[2]    9.48   9.48 1.39 1.28   7.23  11.86 1.02      254      108
##    m1[3]   10.81  10.83 1.15 1.08   8.91  12.76 1.02      272      171
##    m1[4]   12.13  12.14 0.95 0.90  10.54  13.70 1.02      329      306
##    m1[5]   13.46  13.48 0.82 0.78  12.03  14.80 1.01      542      571
##    m1[6]   14.78  14.78 0.81 0.75  13.42  16.12 1.00     1437     1104
##    m1[7]   16.11  16.11 0.92 0.84  14.56  17.62 1.00     2226     1403
##    m1[8]   17.43  17.44 1.11 1.05  15.57  19.22 1.00     1595     1344
##    m1[9]   18.76  18.78 1.34 1.31  16.52  20.91 1.00      933     1132
##    m1[10]  20.09  20.09 1.61 1.60  17.42  22.62 1.00      697      967
##    y1[1]    8.08   8.09 3.98 3.77   1.64  14.43 1.00      762     1178
##    y1[2]    9.53   9.57 3.79 3.74   3.32  15.66 1.00     1078     1292
##    y1[3]   10.74  10.75 3.70 3.68   4.73  16.72 1.00     1506     1788
##    y1[4]   12.12  12.11 3.70 3.53   5.99  18.27 1.00     1851     1877
##    y1[5]   13.43  13.47 3.63 3.45   7.32  19.27 1.00     2014     1694
##    y1[6]   14.81  14.81 3.51 3.42   8.87  20.50 1.00     1988     1673
##    y1[7]   16.15  16.20 3.63 3.53  10.15  22.01 1.00     1775     1835
##    y1[8]   17.61  17.59 3.86 3.74  11.52  23.83 1.00     1967     1866
##    y1[9]   18.74  18.75 3.68 3.45  12.76  24.59 1.00     1828     1674
##    y1[10]  20.05  20.01 3.89 3.80  13.60  26.38 1.00     1448     1646

mdl=cmdstan_model('./ex10.stan') 
fn(mdl,data)
## Initial log joint probability = -110.328 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       21      -10.1859   6.68021e-05   0.000461082           1           1       32    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -10.19
##    b0         5.27
##    b1         1.89
##    s          0.54
##    m0[1]     10.93
##    m0[2]     16.71
##    m0[3]      8.83
##    m0[4]     13.36
##    m0[5]      7.41
##    m0[6]     18.47
##    m0[7]      6.36
##    m0[8]     12.16
##    m0[9]     17.07
##    m0[10]    14.01
##    m0[11]     7.27
##    m0[12]    15.62
##    m0[13]    13.99
##    m0[14]    20.30
##    m0[15]    14.82
##    m0[16]    17.33
##    m0[17]    14.32
##    m0[18]    16.75
##    m0[19]    12.50
##    m0[20]    10.97
##    y0[1]     11.09
##    y0[2]     17.45
##    y0[3]      8.81
##    y0[4]     13.56
##    y0[5]      7.19
##    y0[6]     17.92
##    y0[7]      7.66
##    y0[8]     10.82
##    y0[9]     19.23
##    y0[10]    13.67
##    y0[11]     7.02
##    y0[12]    16.29
##    y0[13]    18.61
##    y0[14]    20.06
##    y0[15]    14.97
##    y0[16]    17.67
##    y0[17]    13.19
##    y0[18]    17.10
##    y0[19]    12.77
##    y0[20]    10.52
##    m1[1]      6.36
##    m1[2]      7.91
##    m1[3]      9.46
##    m1[4]     11.01
##    m1[5]     12.55
##    m1[6]     14.10
##    m1[7]     15.65
##    m1[8]     17.20
##    m1[9]     18.75
##    m1[10]    20.30
##    y1[1]      6.35
##    y1[2]      8.88
##    y1[3]     10.17
##    y1[4]     10.00
##    y1[5]     12.39
##    y1[6]     15.03
##    y1[7]     15.20
##    y1[8]     16.75
##    y1[9]     18.97
##    y1[10]    20.18
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median     sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -12.37 -12.04   1.25 0.99 -14.86 -10.98 1.01      602     1090
##    b0       5.30   5.26   0.50 0.51   4.55   6.20 1.01      474      589
##    b1       1.89   1.89   0.11 0.11   1.71   2.06 1.01      588      702
##    s        0.64   0.62   0.19 0.18   0.38   1.01 1.00      760      837
##    m0[1]   10.97  10.94   0.26 0.24  10.58  11.41 1.00      636      777
##    m0[2]   16.76  16.74   0.30 0.29  16.29  17.26 1.00     2113     1268
##    m0[3]    8.86   8.83   0.34 0.33   8.36   9.45 1.01      501      613
##    m0[4]   13.40  13.39   0.23 0.22  13.06  13.78 1.00     1520     1431
##    m0[5]    7.45   7.41   0.40 0.40   6.85   8.14 1.01      479      568
##    m0[6]   18.52  18.49   0.38 0.36  17.95  19.16 1.00     1510     1094
##    m0[7]    6.39   6.35   0.45 0.45   5.72   7.19 1.01      475      571
##    m0[8]   12.20  12.18   0.23 0.22  11.84  12.61 1.00      903     1011
##    m0[9]   17.12  17.09   0.32 0.30  16.64  17.66 1.00     2023     1251
##    m0[10]  14.05  14.04   0.23 0.22  13.69  14.44 1.00     2139     1435
##    m0[11]   7.30   7.26   0.41 0.41   6.70   8.01 1.01      478      568
##    m0[12]  15.67  15.65   0.27 0.25  15.26  16.11 1.00     2315     1604
##    m0[13]  14.03  14.02   0.23 0.22  13.68  14.42 1.00     2131     1456
##    m0[14]  20.35  20.31   0.47 0.45  19.64  21.16 1.00     1181      928
##    m0[15]  14.86  14.85   0.24 0.23  14.48  15.28 1.00     2326     1578
##    m0[16]  17.38  17.35   0.33 0.31  16.89  17.93 1.00     1950     1204
##    m0[17]  14.36  14.35   0.24 0.23  14.00  14.76 1.00     2232     1458
##    m0[18]  16.80  16.78   0.30 0.29  16.33  17.31 1.00     2105     1268
##    m0[19]  12.54  12.52   0.23 0.22  12.19  12.93 1.00     1044     1169
##    m0[20]  11.00  10.98   0.26 0.24  10.61  11.44 1.00      640      777
##    y0[1]   10.53  10.96  27.11 0.93   7.63  15.48 1.00     1842     1866
##    y0[2]   17.57  16.78  33.14 0.96  12.93  21.03 1.00     1637     1973
##    y0[3]   -2.84   8.86 521.28 1.03   5.45  13.28 1.00     1852     1930
##    y0[4]   14.00  13.43  34.84 0.91   9.73  17.96 1.00     2046     1890
##    y0[5]    7.07   7.41  18.61 1.13   3.34  11.94 1.00     1655     1723
##    y0[6]   17.99  18.53  43.71 1.02  15.18  22.07 1.00     1923     1915
##    y0[7]    6.69   6.35  38.70 1.12   2.14  11.44 1.00     1910     1725
##    y0[8]   10.29  12.14  96.11 0.93   8.12  16.44 1.00     1808     1822
##    y0[9]   20.62  17.08 203.42 0.97  13.22  20.90 1.00     1958     1932
##    y0[10]  14.00  14.04  24.09 0.93  10.41  17.75 1.00     1980     1972
##    y0[11]   4.83   7.27 199.45 1.06   3.03  11.69 1.00     1590     1760
##    y0[12]  15.34  15.67  24.00 1.02  11.86  20.25 1.00     2164     1890
##    y0[13]  15.15  14.02  43.40 0.98  10.06  18.25 1.00     1994     1926
##    y0[14]  20.33  20.35  21.39 1.13  16.15  24.46 1.00     1808     1971
##    y0[15]  14.02  14.84  13.89 0.93  10.22  18.37 1.00     2045     1897
##    y0[16]  17.90  17.36  19.60 1.06  13.24  21.22 1.00     1874     1898
##    y0[17]  14.59  14.32  83.19 0.97  10.07  19.16 1.00     2031     1974
##    y0[18]  16.82  16.81  16.24 1.01  12.77  21.34 1.00     2089     1883
##    y0[19]  12.96  12.52  33.13 0.97   8.49  16.42 1.00     2019     1932
##    y0[20]   6.24  11.03 203.23 0.99   7.42  15.24 1.00     1930     1931
##    m1[1]    6.39   6.35   0.45 0.45   5.72   7.19 1.01      475      571
##    m1[2]    7.94   7.91   0.38 0.37   7.38   8.59 1.01      483      582
##    m1[3]    9.49   9.46   0.31 0.30   9.03  10.04 1.01      523      691
##    m1[4]   11.04  11.02   0.26 0.24  10.66  11.48 1.00      645      777
##    m1[5]   12.60  12.58   0.23 0.22  12.25  12.99 1.00     1074     1197
##    m1[6]   14.15  14.14   0.23 0.22  13.79  14.54 1.00     2167     1458
##    m1[7]   15.70  15.68   0.27 0.25  15.29  16.14 1.00     2311     1604
##    m1[8]   17.25  17.23   0.32 0.31  16.77  17.80 1.00     1983     1251
##    m1[9]   18.80  18.78   0.39 0.37  18.21  19.47 1.00     1440     1107
##    m1[10]  20.35  20.31   0.47 0.45  19.64  21.16 1.00     1181      928
##    y1[1]    6.61   6.39   7.80 1.11   2.41  11.02 1.00     1514     1706
##    y1[2]    7.66   7.90   9.53 0.97   4.15  11.28 1.00     1827     1926
##    y1[3]    9.61   9.51  33.61 0.96   5.58  13.27 1.00     1936     2151
##    y1[4]   10.87  11.05  22.07 0.97   7.59  14.77 1.00     1735     1888
##    y1[5]   12.51  12.57   8.93 0.95   8.44  16.11 1.00     1928     1984
##    y1[6]   14.06  14.13  16.35 0.96  10.91  18.08 1.00     2055     1899
##    y1[7]   15.67  15.68  17.00 0.96  11.35  20.10 1.00     1891     2059
##    y1[8]   17.19  17.23  12.93 1.10  12.43  21.75 1.00     1908     1984
##    y1[9]   19.32  18.82  27.96 1.00  14.71  22.52 1.00     1857     1866
##    y1[10]  21.78  20.35  56.80 1.05  16.13  24.84 1.00     1489     1955

censored data

y i=1-N, y~N(m,s)  
  acutualy        ya i=1-Na
  lower censored  yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
  upper censored  yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)

cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))

P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)

ex11.stan

objective variable is censored

data{
  int N;
  vector[N] x;
  vector[N] y;
  real L;
  int Nl;
  vector[Nl] xl;
  real U;
  int Nu;
  vector[Nu] xu;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
  for(i in 1:Nl)
    target+=normal_lcdf(L | b0+b1*xl[i],s);
  for(i in 1:Nu)
    target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)

xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan') 

mle=mdl$optimize(data=data)
## Initial log joint probability = -6470.48 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       22      -13.5528    0.00361578   0.000438694           1           1       40    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -13.55
##    b0        17.65
##    b1         1.11
##    s          2.73
##    m0[1]     21.39
##    m0[2]     21.89
##    m0[3]     25.42
##    m0[4]     24.32
##    m0[5]     23.65
##    m0[6]     23.82
##    m0[7]     22.15
##    m0[8]     21.65
##    m0[9]     20.24
##    y0[1]     18.06
##    y0[2]     21.57
##    y0[3]     26.36
##    y0[4]     28.30
##    y0[5]     23.82
##    y0[6]     28.16
##    y0[7]     21.27
##    y0[8]     23.55
##    y0[9]     17.85
##    m1[1]     17.92
##    m1[2]     18.99
##    m1[3]     20.07
##    m1[4]     21.14
##    m1[5]     22.22
##    m1[6]     23.29
##    m1[7]     24.37
##    m1[8]     25.44
##    m1[9]     26.52
##    m1[10]    27.59
##    y1[1]     18.16
##    y1[2]     16.36
##    y1[3]     20.64
##    y1[4]     23.39
##    y1[5]     22.29
##    y1[6]     20.96
##    y1[7]     17.38
##    y1[8]     23.67
##    y1[9]     24.31
##    y1[10]    28.30
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -14.41 -14.07 1.48 1.22 -17.74 -12.76 1.01      308      295
##    b0      17.83  17.95 4.59 4.39  10.67  25.19 1.03      227      156
##    b1       1.07   1.04 0.97 0.89  -0.48   2.53 1.03      244      200
##    s        3.89   3.57 1.33 1.06   2.36   6.48 1.00      649      615
##    m0[1]   21.42  21.45 1.74 1.64  18.52  24.16 1.01      290      194
##    m0[2]   21.91  21.96 1.51 1.38  19.42  24.27 1.01      381      402
##    m0[3]   25.31  25.28 2.81 2.44  20.83  29.73 1.02      388      415
##    m0[4]   24.25  24.25 2.02 1.82  21.06  27.47 1.01      562      478
##    m0[5]   23.60  23.61 1.63 1.45  21.05  26.18 1.01      888      710
##    m0[6]   23.77  23.76 1.72 1.53  21.06  26.48 1.01      773      613
##    m0[7]   22.16  22.21 1.42 1.28  19.80  24.44 1.00      477      557
##    m0[8]   21.68  21.72 1.61 1.50  19.03  24.20 1.01      326      299
##    m0[9]   20.32  20.39 2.50 2.42  16.20  24.30 1.02      237      146
##    y0[1]   21.39  21.50 4.63 4.05  13.67  28.75 1.00     1098      901
##    y0[2]   21.65  21.74 4.46 3.95  14.36  28.80 1.00     1567     1703
##    y0[3]   25.23  25.25 4.96 4.27  17.23  32.93 1.01      761      523
##    y0[4]   24.29  24.30 4.57 3.97  16.80  31.31 1.01     1549     1440
##    y0[5]   23.54  23.51 4.60 3.92  16.37  30.68 1.00     2014     1653
##    y0[6]   23.87  23.83 4.50 3.87  16.67  30.99 1.00     1863     1892
##    y0[7]   22.12  22.17 4.32 3.79  15.25  28.81 1.00     1362     1141
##    y0[8]   21.58  21.61 4.43 3.89  14.33  28.53 1.00     1646     1476
##    y0[9]   20.25  20.13 4.76 4.37  12.91  27.89 1.00      806      723
##    m1[1]   18.09  18.19 4.36 4.21  11.30  25.11 1.03      226      156
##    m1[2]   19.12  19.16 3.48 3.34  13.62  24.73 1.03      228      157
##    m1[3]   20.15  20.20 2.63 2.55  15.84  24.32 1.02      235      155
##    m1[4]   21.19  21.25 1.89 1.80  18.06  24.10 1.02      267      126
##    m1[5]   22.22  22.25 1.41 1.26  19.88  24.47 1.00      511      607
##    m1[6]   23.26  23.25 1.48 1.30  20.93  25.59 1.00     1146      872
##    m1[7]   24.29  24.28 2.05 1.84  21.06  27.57 1.01      550      469
##    m1[8]   25.33  25.30 2.82 2.46  20.82  29.80 1.02      386      415
##    m1[9]   26.36  26.32 3.68 3.25  20.47  32.06 1.02      330      347
##    m1[10]  27.40  27.34 4.57 4.07  20.05  34.52 1.02      304      278
##    y1[1]   18.04  18.05 5.92 5.35   8.43  27.41 1.02      382      388
##    y1[2]   19.10  19.15 5.46 4.97  10.15  27.98 1.01      459      433
##    y1[3]   20.01  20.17 4.98 4.28  11.49  27.83 1.01      645      915
##    y1[4]   21.11  21.04 4.63 4.05  13.69  28.35 1.01     1032     1378
##    y1[5]   22.37  22.46 4.29 3.77  15.21  29.10 1.00     1955     1771
##    y1[6]   23.21  23.24 4.33 3.87  16.46  29.97 1.00     2048     1796
##    y1[7]   24.15  24.20 4.52 3.78  16.81  31.51 1.00     1846     1274
##    y1[8]   25.38  25.49 4.98 4.47  17.19  33.30 1.01     1012      875
##    y1[9]   26.42  26.59 5.64 4.89  17.31  35.16 1.01      660      705
##    y1[10]  27.54  27.41 6.22 5.35  17.77  37.71 1.01      494      787

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

data=list(N=n,x=xya$x,y=xya$y,
          L=L,Nl=nl,xl=xyl$x,
          U=U,Nu=nu,xu=xyu$x,
          N1=n1,x1=x1)
mdl=cmdstan_model('./ex11.stan') 

mle=mdl$optimize(data=data)
## Initial log joint probability = -68.8876 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       25      -18.2624    0.00137193   0.000255376           1           1       30    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -18.26
##    b0         9.66
##    b1         2.90
##    s          3.32
##    m0[1]     19.47
##    m0[2]     20.79
##    m0[3]     30.04
##    m0[4]     27.17
##    m0[5]     25.39
##    m0[6]     25.85
##    m0[7]     21.48
##    m0[8]     20.15
##    m0[9]     16.46
##    y0[1]     17.57
##    y0[2]     17.99
##    y0[3]     28.02
##    y0[4]     29.44
##    y0[5]     23.86
##    y0[6]     29.86
##    y0[7]     22.77
##    y0[8]     22.97
##    y0[9]     18.34
##    m1[1]     10.37
##    m1[2]     13.19
##    m1[3]     16.01
##    m1[4]     18.83
##    m1[5]     21.64
##    m1[6]     24.46
##    m1[7]     27.28
##    m1[8]     30.10
##    m1[9]     32.92
##    m1[10]    35.74
##    y1[1]     10.95
##    y1[2]     10.25
##    y1[3]     16.16
##    y1[4]     22.28
##    y1[5]     22.72
##    y1[6]     24.66
##    y1[7]     32.29
##    y1[8]     29.59
##    y1[9]     29.18
##    y1[10]    33.12
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,'m',ch=T)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -18.82 -18.39 1.52 1.20 -21.89 -17.25 1.00      355      466
##    b0       6.89   7.75 4.47 3.50  -1.36  12.12 1.01      303      306
##    b1       3.52   3.33 0.91 0.69   2.47   5.14 1.01      317      358
##    s        4.77   4.38 1.84 1.28   2.84   7.99 1.00      512      567
##    m0[1]   18.79  19.01 1.87 1.50  15.45  21.37 1.01      417      472
##    m0[2]   20.40  20.52 1.64 1.31  17.57  22.76 1.01      519      592
##    m0[3]   31.63  31.19 2.60 2.03  28.52  36.15 1.00      589      664
##    m0[4]   28.15  27.87 1.93 1.57  25.75  31.44 1.00      961      778
##    m0[5]   25.99  25.81 1.62 1.36  23.83  28.71 1.00     1433      945
##    m0[6]   26.54  26.33 1.69 1.42  24.34  29.37 1.00     1321      938
##    m0[7]   21.23  21.31 1.56 1.27  18.59  23.48 1.01      610      647
##    m0[8]   19.62  19.78 1.74 1.42  16.58  22.07 1.01      459      531
##    m0[9]   15.14  15.56 2.56 2.08  10.49  18.36 1.01      334      360
##    y0[1]   18.82  18.92 5.44 4.60   9.47  27.51 1.00     1672     1002
##    y0[2]   20.36  20.37 5.31 4.66  11.87  28.97 1.00     1914     1464
##    y0[3]   31.80  31.34 5.93 5.15  23.39  41.99 1.00     1508     1213
##    y0[4]   28.02  27.68 5.45 4.52  19.94  36.66 1.00     1881     1875
##    y0[5]   25.87  25.83 5.26 4.49  17.84  34.12 1.00     2092     1610
##    y0[6]   26.34  26.17 5.28 4.62  18.19  34.91 1.00     1763     1500
##    y0[7]   21.37  21.30 5.44 4.56  13.03  29.60 1.00     1846     1633
##    y0[8]   19.79  20.04 5.19 4.55  11.64  27.80 1.00     1628     1157
##    y0[9]   15.13  15.60 5.53 4.94   4.95  23.11 1.01      980      983
##    m1[1]    7.75   8.58 4.26 3.32   0.04  12.78 1.01      304      307
##    m1[2]   11.17  11.82 3.45 2.68   4.92  15.31 1.01      310      328
##    m1[3]   14.59  15.01 2.67 2.17   9.68  17.92 1.01      328      364
##    m1[4]   18.01  18.25 2.00 1.64  14.39  20.73 1.01      388      417
##    m1[5]   21.44  21.51 1.54 1.25  18.84  23.68 1.01      639      657
##    m1[6]   24.86  24.73 1.52 1.29  22.71  27.36 1.00     1439      867
##    m1[7]   28.28  27.99 1.95 1.58  25.87  31.61 1.00      937      730
##    m1[8]   31.70  31.25 2.62 2.04  28.57  36.23 1.00      585      664
##    m1[9]   35.12  34.51 3.38 2.63  31.10  41.11 1.01      466      488
##    m1[10]  38.54  37.72 4.19 3.27  33.58  45.84 1.01      415      472
##    y1[1]    7.74   8.32 6.43 5.59  -3.67  16.94 1.00      799      650
##    y1[2]   11.05  11.65 6.17 5.23  -0.02  19.76 1.00      690      825
##    y1[3]   14.50  14.89 5.72 4.97   4.60  23.24 1.00      968     1022
##    y1[4]   17.96  18.26 5.66 4.69   8.44  26.32 1.00     1226     1235
##    y1[5]   21.56  21.59 5.32 4.51  13.26  29.97 1.00     1332     1572
##    y1[6]   24.98  24.83 5.35 4.37  16.82  34.16 1.00     1656     1233
##    y1[7]   28.15  27.85 5.43 4.85  20.00  37.13 1.00     1814     1326
##    y1[8]   31.70  31.19 5.74 4.69  23.53  42.43 1.00     1351     1023
##    y1[9]   35.00  34.47 6.31 5.31  26.37  45.51 1.00      901      932
##    y1[10]  38.47  37.82 6.32 5.38  29.82  49.32 1.00      885      915

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

appendix: bimodal distribution, mixed normal distribution

stan

data {
  int<lower=0> N;
  real y[N];
}
parameters {
  real<lower=0, upper=1> theta; //ratio of mix
  real mu1;
  real mu2;
  real<lower=0> sigma1;
  real<lower=0> sigma2;
}
model {
  for (n in 1:N) {
    target += log_mix(theta,
                      normal_lpdf(y[n] | mu1, sigma1),
                      normal_lpdf(y[n] | mu2, sigma2));
  }
}

EM algorithm

y0=rnorm(100,0,1)
y1=rnorm(100,-5,1)
y2=rnorm(100,5,1)
y=sample(c(y0,y1,y2),100)
density(y) |> plot()

library(mclust)

rst=Mclust(y)
summary(rst)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust E (univariate, equal variance) model with 3 components: 
## 
##  log-likelihood   n df  BIC  ICL
##            -242 100  6 -511 -512
## 
## Clustering table:
##  1  2  3 
## 35 34 31
rst$parameters
## $pro
## [1] 0.356 0.334 0.310
## 
## $mean
##       1       2       3 
## -4.8267  0.0605  5.5018 
## 
## $variance
## $variance$modelName
## [1] "E"
## 
## $variance$d
## [1] 1
## 
## $variance$G
## [1] 3
## 
## $variance$sigmasq
## [1] 0.841
plot(rst)